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ABSTRACT 



A simple primitive-equation model for numerical 
prediction of wind fields and temperature fields at BOO and 
400 mb is developed and programmed. Niimerical experiments 
are conducted with this model in an effort to suppress the 
spurious instability encountered in solution of the 
primitive equations. The model proves to be very sensitive 
to spatial finite-differencing schemes. The schemes 
developed in this study tend to suppress the instability 
resulting from truncation and round-off error. Instability 
resulting from unrealistic boundary conditions for a 
limited forecast region is reduced by not allowing 
advection across the boundaries and appears to be con- 
trolled by the finite-differencing schemes used in this 
study. 
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Professor George J. Haltiner of the U. S. Naval 
Postgraduate School for his guidance, contributions and 
encouragement in this work. Appreciation is also expressed 
to the personnel of the U, S. Navy Fleet Numerical Weather 
Facility for their cooperation during the preparation of 
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Introduction. 



It is generally recognized that solution of the 
primitive -equations offers the most general and logical 
approach to numerical prediction. Previous attempts at 
solution of the primitive equations by numerical methods 
have met with limited success. However, interest has 
recently been stimulated in the primitive equations by 
unpublished reports indicating that Arakawa [yj has 
succeeded in developing a stable solution. The purpose 
of this study is to develop and program a simple 3- level 
primitive equation model and perform some numerical 
experiments. The main effort in these experiments is 
devoted to the suppression of spurious instability 
encountered in numerical solution of the primitive 
equations. It is hoped that this instability can be 
suppressed by a finite-differencing scheme that will in 
the mean conserve the sum of potential and kinetic energies. 
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The Model 



Phillips [ll3 describes a coordinate system where the 
vertical coordinate p in the x, y, p, t-system is replaced 
by an independent variable V, where ^ = 
relation then holds, where Q can be x, y, or t: 

/J-] = i ilL . 

U&/p -nr (II 

The subscript p denotes derivatives taken along a constant 
pressure surface. The variable ST ranges monotonically from 
zero at the top of the atmosphere to unity at the ground. 

If the above coordinate system is utilized the 
horizontal equation of motion becomes 

(^V +• V(|) V7T + 

.^t aT TTdT 

O. (21 



In this study it is assumed that the top of the 
atmosphere is at 200 mb. The atmosphere is divided 
vertically by five levels as follows: 



II 

o 


p=200 mb 


^ =l/4_._ 


p=400 mb 


=1/2 _ . _ . 


p=600 mb 


=3/4 


p=S00 mb 


^=1 


SFC 



-Pp. The following 



Vi' 



- 2 - 




- « • • • 





I • ^ 



[ 



For this model and will be approximated as follows 






(3) 



X 

Vertical velocities are identically zero at the top and 

♦ • 0 

bottom of the atmosphere (^=0,^=0), reducing ^ and 

* ' 

• ... 

2- 



If the preceding simplifications are utilized and 
equation (2) is applied at the and '^levels, the 
following equations ensue 



iK ^ V‘p, 

dt SAT 



VT V/7-t If t fkxX ~ O 

rr aT 

+ ('^•'=7)'^ + Vx (V 3 -V ) + - 

^ u a2iV” 

S V7T+ F t f = O . 

V- ar 



It is assumed in (4) and (5) that 

^Vi n. . 

'T" d ^ r:r 
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Moreover the term ^ ^ in (4) may be replaced by 




0 



Since ^=0 and =1, 

= <kSL - 4>1 J 

aT 

and equation (4) becomes 

+^-'7h+i^(Vi-^) + v4), t(4>, - 

a u oLoy 



if. i- KX£'^~ O- 

Similiarly at the level 

^ ~ ^ ~ 03 

av- 

and equation (5) becomes 

|^4^.v)'^+^(V3~'v) -t- V4^t [d>i- 

+ 1<X£'^=: O. (7) 



Equations (6) and (7) may also be expressed in component 
form as follows: 
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U.,^1 + - Xl S ]-SjX3“X) 

_ < 3 / rfX J 

C“'j|'+<§+ ’i'’'']-^v<’'j-<' -if 



(g) 



<)t 



(9) 



^3 - - K -'''3 % 1 - ^ (^3' 






( 10 ) 



i^ = ^ Vj s^yj + 7W7».3 l-^iL fx ~'^)-sMl3 - 

L <5^ J 2^7- 

[^3 )]:j^ - 5^ (II 

The equation of continuity expressed in ^ coordinates 
becomes 



V • frrv) + i-iiTT) - O 

c>t 



( 12 ) 



This equation can be combined with the thermodynamic 
equation 



a 6 -u V- ve-t- ~ o 

s 



(13) 
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to yield 



ci(ll Q) j- V {TTQ'y) + 'T-) c: o. <!'*> 

J t 

Taking the dot product between fl V =ind equation (2) 
leads to 

jr JTV- vv"'4 

^ ^~t X ^ 

77 V‘ '7(p - TV* 777 - o# (15) 

Multiplying (12) by (p and utilizing the relationship 



Try-v<ip = y-(Trvc)>)- (jfv-dT'y} 



results in 



TfV- V 4 > = V- (ttV 4 >) +■ -I- 

c) 77 '"7 



Substituting the identities 



(-i» 



■AT 



and 



into (15) and adding (15) to (12) yields a kinetic energy 
equation of the form 
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i_( X v"1+ V- (■? v+ 'y ) -f -^ ( -XT v"" X i- ]t<Pt) -\ 



i_ (4>T)ilT _(Xinr -/-XV-V7T-/"77'7")i^ - O. (16) 

c)^ r^t <^'t 



The first law of thermodynamics applied to adiabatic 
conditions is 



Cp dr — c<. 4£ - o 
<AC <st 



Substituting p= TTV" ‘•’Pt expanding the total derivatives 
into Eulerian form for coordinates leads to 

i_ (CpT) + V- VfCpT) + - cK r ^JLSE) f- 

dt Ldt 

\y. vfrr v) + at Mils:) 1 = o ■ 

d'T- J 

Multiplying this equation by yp and equation (12) by 
CpT, then adding the resulting equations gives the 
potential energy equation 



^(TCpl) + V-(rrCpTVj+ AiFCpTT)-7rd,l mi3)4. 

=>t 

V-V(7t^)+^4IjI3^1 -O . (17) 

J 
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Combining equations (16) and (17) gives 

i_(TTCVT+ n; \/' V-f/TCp IVt jrv' Vt 7T<i) 'y) -t- 
c)t ^ 

J_(TTCpT'HlE V*^ + 77 c( 3V-)+ i]T4j:iS:') = cP. msi 

a- 

Thompson outlines the following three main aspects 

to be considered in solving the primitive equations : 
formulation of boundary conditions, formulation of initial 
conditions and numerical integration of finite difference 
equations. 

The method of numerical integration utilized in this 
study will be the method of repeated time extrapolation 
outlined by Thompson and shown to be computationally 
stable for a simple linear model with a grid interval the 
size of the one used in this study if the time interval 
is about 10 min. or less. This method employs a centered 
difference of the form 




Initial conditions for solution of the primitive 
equations present a real problem. For meteorologically 
significant large-scale motion there is a near balance 
between the coriolis and pressure forces. This means the 
resultant acceleration is an order of magnitude smaller 
than either of these two forces. Small errors in the 





I f 












i i 



^ 1 ^ 



* ll«^ 






measurement of initial wind and pressure fields may lead 
to intense divergence fields which, in turn, give large 
pressure changes, thus creating rapid accelerations., 

These spurious accelerations create fictitious gravity 
waves. To exclude these waves initially in this study, 
a form of non-divergent initial winds is used. However, 
this does not necessarily exclude gravity waves for all 
time . 

The forecast region used in this study covers the 
major portion of the Northern Hemisphere. Boundary 
conditions for a limited region such as this are difficult 
to define since the proper type of boundary conditions 
depends on the mathematical character of the primitive 
equations, which have a shifting dual nature, sometimes 
elliptic and sometimes hyperbolic. In this paper all 
parameters are assumed constant with respect to time on 
the boundaries. It is recognized that this assumption 
may lead to the existence of spurious development near the 
boundaries; however, it is felt that if the spatial finite- 
difference schemes are such that in the mean potential and 
kinetic energies are conserved, any boundary instability 
will tend to be averaged out over the forecast region. This 
averaging process will also tend to suppress spurious gravity 
waves resulting from truncation and round-off errors. 

It has been shown by application (see, for example, 
Haltiner [6 J ) of Gauss’ equation that when an equation 



-9 



of the same type as (IS) is integrated over a closed 
volume the local time change of the sum of kinetic and 
potential energy vanishes. It follows that any finite- 
difference scheme should in the mean conserve the sum of 
the potential and kinetic energies. Specifically, this 
consistency should be achieved separately with respect to 
both vertical and horizontal differencing. 

The vertical consistency is developed as follows. 
First equation (15) is applied at the level. The 
vertical differencing is then performed, leading to 



is approximated by a one-sided difference, A similar 
application of equation (1?) results in 






(19) 



since ^ =0, ^ =0 and derivative 




J_(7TC,T7) + V- (VCpT, t ^(VTTCpT.') 

/i\r 








( 20 ) 
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Adding equations (19) and (20) and comparing the result 
with equation (IS) leads to the conclusion that if the 
consistency discussed in the previous paragraph is to be 
maintained, certain terms of equation (19) must cancel 
corresponding terms of equation (20)o This inference 
leads to the following equivalences : 

7T5L_ r? T; ^0, - 0a = Cp C - 
TTV+f^ 




A similar procedure at the level utilizing the 
identity 




leads to 



TJ V. R Tj ^ Cpa + 05 - a. 04 
rr^+Py 

and 

4-03-Cp T^) - Q. 

From the previous equations the following relationships can 
be obtained. 




Substituting P =200 mb, p.=400 mb, p =S00 mb 
t 3 

k=0.2S6j/g°k, and C =1.003 j/g°k in the previous system 

r' 
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of equations leads to 



?2= .643 T]_ + .410 T 3 
4^1 = $^H+ .442 RTi + .530 RT 3 
0a = “. 05 S RTi + .530 RT 3 
03 = 04 + . 05 s RTi + ,22 RT 3 



(21) 

( 22 ) 

(23) 

(24) 



If these relationships between geopotential and temperature 
are maintained then the vertical differencing scheme will 
not violate the energy conservation equations. 

With respect to simple derivatives of the form 



involved in (15), the centered difference approximation 



will obviously vanish when summed over interior points of 
the grid and thus would in the mean conserve potential and 
kinetic energy. However this does not take care of the 
boundaries, nor terms which are not of the simple form 
shown above. 

Since numerical integration is to be performed by 
repeated time extrapolation, the velocity tendencies (S) - 



have been extrapolated forward in time. This will require 







(11) will have to be computed from 




tendency equations for JT ~ and 

in terms of T through equations (22) - (24). 

Integration of (12) leads to the approximate form 





(25) 



^12- 



From (12) 



US. - V Cttv) 

at- 



substituting for.c^JQT from ( 25 ) gives 

^[s/-(rrv,~TTV,)l . 



The relationship 

-TTci, =: M 



Utilized in equation (I?) yields for the VJ" level 

TTCp^i- Cpsl air-+ v-('v\^)/tCpTr'i/-\7r + 

J b J 



00 - <t>/ 1 ^4- V- yn-7 f vr - 4^/ o. 

av Ut j 



( 26 ) 



The bracketed terms of equation (26) may be replaced by the 
following relationship 

^ -t- V-fTTO/; = -21S-. 



Dividing (26) by 7T Cp gives 

UL + V-VT; =: O . 



(27) 









A similar approach at the level will result in the 
following tendency equation 



VT3 + . t2 - .033T1 









=0 



{2B) 



A3' 



based on equations (21) - (24) with <|)<y assumed zero. 

These tendency equations complete the system required 
for repeated time extrapolation. 
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3. Numerical experiments and results. 

1. Energy Considerations 

In the basic model, the energy balance with respect 
to vertical differencing is maintained through the 
relationship between temperature and geopotential developed 
in Section 2. With respect to horizontal differencing, 
the problem may be considered in two parts, 

A. internal energy balance, 

B. boundary effects. 

For the former problem, partial derivatives of a 
"closed" form, such as^_£i^, can be maintained in balance 
through use of a simple central differencing scheme. As 
the computation proceeds across the grid, each term is 
added and subtracted once at each grid point, with no 
change in total energy. However, a few derivatives with 
variable coefficients must be evaluated in the horizontal 
acceleration equations, (S) through (11), and also in the 
temperature-change and stability^change equations. In 
these terms, of the form v_^, the coefficient varies as 
the grid is traversed, so that the weighting of the term 
at a specific grid point changes between addition and sub- 
traction in central differencing. This can lead to a change 
in total energy across the field, with resultant build- 
up of instability. 

The following approaches were tested on these "open" 
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forms : 



a» Terms left in this form» 

This lead to a comparatively short set of forms « For 
example, the bracketed portion of equation (S) was 
simplified to ^1'^ ^1 =■ ^1^ ^1 •“ v-.fo However, this form did 
not fulfill the requirements for energy balance, and the 
model became unstable in a very short time with overflow 
occuring in less than four hours » 

b» Terms placed in a ’’closed” form as a sum of squares « 
The bracketed form of equation {S} now become 

This form was a distinct improvement on the ’’open” form; 
however, variable coefficients involved in7|j, were not 
affected. In any case, there was some instability and this 
model overflowed in about six hours. 

Co Open form with multipliers replaced by average 
values o 

In this model, the equations of the first approach 
were used, but a "central average” replaced the variable 
coefficient in each case. For example, Vg)u became, in 
finite difference notation: 

(Vf.T-tA +Vi(T-h ■) (U'C>T-tO- I 

2 2^y 

This model gave the best results of the three, and was 
used for the remaining experiments. 

The problem of boundary effects proved much more 
difficult. Due to computer storage limitations, this model 
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was programmed using an octagonal grid whose southern 
boundaries are as much as from the equator c. Advection 
at these boundaries leads to major problems » 

The following boundary forms were tested: 

a. All time derivatives set equal to zero along the 
boundaries. Since the boundary values were held constant, 
while values at adjacent points were allowed to vary, 
large gradients resulted, giving large-amplitude gravity 
waves which moved inward from all sides at high speed. 

(See Fig. 1) 

b. Boundary values averaged at short time intervals. 

Several smoothing techniques were tried which 

involved only the outside three grid points. Since a 
gravity wave cannot cover the distance between two grid 
points in less than one hour of forecast time, this 
averaging was carried out once each hour. The most 
successful technique involved placing an average of the 
three outer points at the next-to-edge point, then averaging 
this value with that at the outside point and placing the 
result at the outside point. For example; 

(Boundary) A B C 

becomes A+(A+B+C)/3 (A+B+C)/3 C 

2 

This method restricted the instability due to boundary 
effects to the outer two grid points; however, in spite of 
the heavy averaging, large differences developed between 
the outer three grid points after about five hours fore- 
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cast time. 



c. All horizontal velocities set equal to zero along 
the boundaries. 

Since the instability developed along the boundaries 
can be traced to energy advected into or out of the forecast 
area, zeroing the horizontal winds at the boundaries, 
thereby zeroing the advection, should solve the problem. 

In practice, it was found that zeroing the horizontal wind 
components at both forecast levels at the outer two grid 

I 

points was effective in controlling boundary-produced gravity 
waves. Surprisingly, this did not lead to instability in 
the wind-component fields themselves in this model. This 
method of boundary condition control was used throughout 
the remaining experiments. 











% 



2. Friction Terms 

Two friction terms were developed for use with the 

horizontal acceleration equations (^) through (11). For 

horizontal friction, a simple Laplacian multiplied by lO^^ny^ec- 

was subtracted at each time step. For vertical friction, 

the terms developed by Smagorinsky 5j were used, with an 

average skin friction coefficient of 2.2 x 10~3 after 

Cressman [3]- According to Smagorinsky, the vertical 

friction at any level may be computed from 

vF = -go[3^= g_ n^+i "X-i)* 

^ cy p ^p 

Assuming 7o=0, 

; vF=g_(%-7i) 

■L /I p J AP 

Vertical stress, y' , is taken as a function of the vertical 
wind shear, so that 

or since K is proportional to the magnitude of the wind, 

1^1 4 f- is a constant for any roughness length 

=(K iviv| 

^d depends on roughness length, and from Cressman 's paper 
2 -3 

an average value of 2.2 x 10 was used. Level 4 was taken 
to be the top of the friction level for this tern, with 
I vj^^ estimated as a fraction of the vectorially extrapolated 
wind by the equation 

ivlf = 0.36 ( (1.192u--0.192u )^+ (1.192v.-0.192vJ^) . 

* * 4 J- 
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For\y^ the cross isobar angle at level 4 was taken to be 
a function of latitude alone by the equation 
ctncT= l“ 10 ^l 2 f j f in 

This gives a cross-isobar angle ranging from near zero 
at the pole to approximately 35° at the base of the grid 
(ISN), which is in close agreement with average values* 

The magnitude of was obtained from the geostrophic 
approximation and the A field j which in turn was estimated 
by a simple extrapolation from (p^ and (pj j namely 

<P-f =1,192 - 0 , 192 , 

For Tv the equation, 

X- -(n)2 kA (^"Vo) 

^ ^p 1 i 

can be readily obtained from the basic T* equation. Here 
2 , where z=5»236 km, the standard depth of the 
SOO- to 400-mb layer. For (^K) 2 » the Rossby and 
Montgomery estimate of 500 gm/m sec was used. 

This series of equations was programmed and tested 
on the model. However, the terms involved were so small 
compared to the rest of the model, and the computations 
so time-consuming (approximately 50^ increase in computing 
time), that time was not available for a complete evaluation. 
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3. Prognosis of Surface Pressure 



The prognosis of surface pressure was one of the most 
difficult problems encountered because the surface pressure 
change of the basic model was apparently very sensitive to 
small errors in velocity aloft. In fact, the surface 
pressure field served as a diagnostic chart, since any 
instability was apparent there long before it appeared 
aloft. 

Four different methods were tried for obtaining the 
prognostic surface pressure field; 

a. Expansion of the basic equation (25). 

The resulting pressure-change field was extremely 
sensitive to changes in velocity divergence at levels 1 and 
3, and large pressure changes occurred in two to three hours. 

b. Averaging around the grid point. 

In this approach, the products of the pressure value 
and the sum of the wind components at levels one and three 
were differenced across the point being computed. In 
finite difference form, 



This model was not quite as sensitive to divergence 
errors, but still gave unreasonable surface pressure changes 
over periods as short as 6 hours. 



c. Computation of surface pressure directly from the 
fields. 
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In this case, the deviation of the surface pressure 
from standard was estimated at each time step from the 
deviation of the (j) fields from standard at SOO and 400 mb 
and the hydrostatic equation^ i.e., if is defined as 

/ / 

and if we let ■=. > ^^en from the hydrostatic 

equation, and p=Pgtd P- 

This model was programmed and computed. However, 
since the (p fields are based on temperature and the basic 
inputs to this model were the temperature patterns at 400 
and SOO mb, the resultant surface pressure field lacked 
some of the usual detail. 
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Conclusions and suggestions for further study 

In this series of experiments j three finite difference 
forms, three boundary forms, and three surface pressure 
forms were tested. It was found that this model is 
extremely sensitive to finite difference forms used and 
boundary conditions assumed. However, when a new form gave 
a closer approximation to energy balance, stability was 
improved in every case, with the improvement noted being 
in direct proportion to the approach to energy balance. 

It is the opinion of the writers that a stable 
primitive-equation model of this type can be achieved, and 
will eventually become the standard form of numerical 
prognosis. However, faster computers with much larger 
storage capacity must be used to realize this goal. At 
least two more levels of temperature input are needed to 
delineate the geopotential fields clearly, and the grid 
should be expanded to cover at least a whole hemisphere 
to reduce boundary problems. Even with high speed magnetic 
drum storage , such a model would have a prohibitively long 
running time on a computer such as the CDC 1604; but it 
might be completely feasible on the next generation of 
computers . 

For the present, further experiments with surface 
pressure forms are suggested, possibly using some simpli- 
fied version of the equations currently under development 
by Arakawa • Also, the friction terms could be 
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streamlined and tested, and a surface terrain term included. 
However, any major revision of this model will require the 
use of a larger, faster computer than any currently 
available . 
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LIST OF FIGURES 

Figure 

1. Section of 6-hour SOO-mb temperature prognosis with 
no boundary correction. 

2. Section of 12-hour SOO-mb temperature prognosis, 
edge-averaged. 

3. Section of iS-hour SOO-mb temperature prognosis, 
boundary advection zeroed. 
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